rm(list=ls())
library(tidyverse)
library(reshape2)
library(e1071)

Support vector machines

  1. Suppose a categorical label: \(y\in\{0,1\}\) and a data generating process: \(x|y=0\sim\mathcal{N}(0,1)\) and \(x|y=1\sim\mathcal{N}(2, 1)\). Find the optimal classification boundary. In doing so, determine the Bayes optimal classification accuracy. In both cases, assume categories are equally as likely to occur in the data.

Drawing the two densities

n <- 100
x <- seq(-2, 5, length.out = n)
y0 <- map_dbl(x, ~dnorm(.))
y1 <- map_dbl(x, ~dnorm(., 2, 1))
tibble(x, y0, y1) %>% 
  melt(id.vars="x") %>% 
  ggplot(aes(x, y=value, colour=as.factor(variable))) +
  geom_line()

We see that the the densities crossover at a point. Let’s find what it is.

\(\mathcal{N}(x|0,1) = \mathcal{N}(x|2,1)\)

which, by symmetry, occurs when \(x=1\).

So the optimal boundary is predict \(y=0\) when \(x < 1\) and \(y=1\) when \(x>1\). At the boundary, predict either with probability \(1/2\).

The optimal classification accuracy is \(Pr(x<1|y=0) = \int_{-\infty}^{1} \mathcal{N}(x|0,1) \approx 84\%\).

  1. Assume that \((x_1,x_2)'|y=0\sim\mathcal{N}((0,0)', I)\) and \((x_1,x_2)'|y=1\sim\mathcal{N}((1, 1)',I)\). Simulate 100 draws from each category. Use the parametric form of each conditional density to classify each point and visualise these; in doing so, calculate the optimal classification rate on these data.
library(mvtnorm)
rmvrnorm2D <- function(n, mux, muy, sigmax, sigmay, rho){
  return(rmvnorm(n, c(mux, muy),
                 matrix(c(sigmax^2, sigmax * sigmay * rho,
                          sigmax * sigmay * rho, sigmay^2),
                        ncol = 2)))
}

n <- 100
x0 <- rmvrnorm2D(n, 0, 0, 1, 1, 0) %>% 
  as.data.frame() %>% 
  mutate(class=0)
x1 <- rmvrnorm2D(n, 1, 1, 1, 1, 0) %>% 
  as.data.frame() %>% 
  mutate(class=1)

# true classes
x <- x0 %>% 
  bind_rows(x1)
x %>% 
  ggplot(aes(x=V1, y=V2, colour=as.factor(class))) +
  geom_point()

Calculate the probability density for each point. In doing so, this yields a probability of belonging to each class (where we know each are equally likely):

\(Pr(y=1|x_1,x_2) = \frac{p(x_1,x_2|y=1) Pr(y=1)}{p(x_1,x_2)}=\frac{p(x_1,x_2|y=1) Pr(y=1)}{p(x_1,x_2|y=1) Pr(y=1) + p(x_1,x_2|y=0) Pr(y=0)}\)

dmvrnorm2D <- function(x, mux, muy, sigmax, sigmay, rho, log=FALSE){
  return(dmvnorm(x, c(mux, muy),
                 matrix(c(sigmax^2, sigmax * sigmay * rho,
                          sigmax * sigmay * rho, sigmay^2),
                         ncol = 2),
                 log))
}

p0 <- map_dbl(seq_along(x$V1), ~dmvrnorm2D(x[., 1:2], 0, 0, 1, 1, 0))
p1 <- map_dbl(seq_along(x$V1), ~dmvrnorm2D(x[., 1:2], 1, 1, 1, 1, 0))

df <- x %>% 
  mutate(p0=p0, p1=p1) %>% 
  mutate(class_pred=if_else(p1 > p0, 1, 0)) %>% 
  mutate(correct=if_else(class==class_pred, 1, 0))

mean(df$correct)
[1] 0.78

Visualise incorrect predictions

df %>% 
  ggplot(aes(x=V1, y=V2, colour=as.factor(correct))) +
  geom_point()

  1. What is the optimal classification boundary here?

The optimal classification boundary is given by:

\[\begin{equation} \mathcal{N}(x_1,x_2|(0,0)', I) = \mathcal{N}(x_1,x_2|(1,1)', I) \end{equation}\]

Solving these equations, we find that: \(x_2= 1 - x_1\). Plotting this on the data distribution:

x %>% 
  ggplot(aes(x=V1, y=V2, colour=as.factor(class))) +
  geom_point() +
  stat_function(fun = function(V1) 1 - V1, colour="black") +
  coord_fixed() +
  scale_y_continuous(limits = c(-4, 4)) +
  scale_x_continuous(limits = c(-4, 4))

  1. Apply (linear) SVM to the dataset you generated in the previous question. What’s its predictive accuracy on your dataset? How close does it get to the optimal classification boundary?

Not bad but overfits a bit.

df <- df %>% 
  mutate(class=as.factor(class))
svmfit <- svm(class ~ ., data = df %>% select(V1, V2, class), kernel = "linear", cost = 10, scale = FALSE,
              probability=TRUE)

## predictive accuracy
class1 <- predict(svmfit, df)
df %>% 
  select(class) %>% 
  mutate(class_pred=class1) %>% 
  mutate(correct=if_else(class==class_pred, 1, 0)) %>% 
  summarise(mean(correct))

# plot classification boundary
test_df <- expand_grid(V1=seq(-4, 4, length.out = 100), V2=seq(-4, 4, length.out = 100))
classes <- predict(svmfit, test_df)
test_df %>% 
  mutate(class_pred=classes) %>% 
  ggplot(aes(x=V1, y=V2, colour=as.factor(class_pred))) +
  geom_point() +
  stat_function(fun = function(V1) 1 - V1, colour="black") +
  coord_fixed() +
  scale_y_continuous(limits = c(-4, 4)) +
  scale_x_continuous(limits = c(-4, 4))

How do predicted class probabilities compare with true ones?

class_probs <- predict(svmfit, df, probability = T)
  1. Fit an SVM with a radial basis function kernel. What’s its predictive accuracy? What do its classification boundaries look like?
df <- df %>% 
  mutate(class=as.factor(class))
svmfit_radial <- svm(class ~ ., data = df %>% select(V1, V2, class), kernel = "radial", cost = 10, scale = FALSE)

## predictive accuracy
class1 <- predict(svmfit_radial, df)
df %>% 
  select(class) %>% 
  mutate(class_pred=class1) %>% 
  mutate(correct=if_else(class==class_pred, 1, 0)) %>% 
  summarise(mean(correct))

# plot classification boundary
test_df <- expand_grid(V1=seq(-4, 4, length.out = 100), V2=seq(-4, 4, length.out = 100))
classes <- predict(svmfit_radial, test_df)
test_df %>% 
  mutate(class_pred=classes) %>% 
  ggplot(aes(x=V1, y=V2, colour=as.factor(class_pred))) +
  geom_point() +
  stat_function(fun = function(V1) 1 - V1, colour="black") +
  coord_fixed() +
  scale_y_continuous(limits = c(-4, 4)) +
  scale_x_continuous(limits = c(-4, 4))

  1. Compare the predictive accuracy of each SVM model on a new dataset drawn from the same distribution. Which SVM performs best?

It should be the linear one as the radial one overfits the data.

n <- 100
x0 <- rmvrnorm2D(n, 0, 0, 1, 1, 0) %>% 
  as.data.frame() %>% 
  mutate(class=0)
x1 <- rmvrnorm2D(n, 1, 1, 1, 1, 0) %>% 
  as.data.frame() %>% 
  mutate(class=1)
x <- x0 %>% 
  bind_rows(x1)
class_linear <- predict(svmfit, x)
class_radial <- predict(svmfit_radial, x)

x %>% 
  mutate(class_linear=class_linear,
         class_radial=class_radial) %>% 
  summarise(mean(class==class_linear), mean(class==class_radial))
  1. Generate data for one class by sampling \((x_1,x_2)'|y=0\sim\mathcal{N}((0,0)', I)\). For the other class sample from \((x_1,x_2)'|y=1\sim\mathcal{N}((0,0)', 2 I)\) such that \(x_1^2 + x_2^2 \geq 3^2\).
n <- 100
x0 <- rmvrnorm2D(n, 0, 0, 1, 1, 0) %>% 
  as.data.frame() %>% 
  mutate(class=0)

f_sample_constrained <- function(r) {
  x_test <- rmvrnorm2D(1, 0, 0, 2, 2, 0)[1, ]
  while((x_test[1]^2 + x_test[2]^2) < r^2)
    x_test <- rmvrnorm2D(1, 0, 0, 2, 2, 0)[1, ]
  return(x_test)
}

x1 <- matrix(nrow = n, ncol = 2)
for(i in 1:n)
  x1[i, ] <- f_sample_constrained(3)

x1 <- x1 %>% 
  as.data.frame() %>% 
  mutate(class=1)

x <- x0 %>% 
  bind_rows(x1) %>% 
  mutate(class=as.factor(class))
x %>% 
  ggplot(aes(x=V1, y=V2, colour=class)) +
  geom_point()

  1. Fit a classifier using radial basis functions and visualise the decision boundary
svmfit_radial <- svm(class ~ ., data = x, kernel = "radial", cost = 10, scale = FALSE)

test_df <- expand_grid(V1=seq(-6, 6, length.out = 100), V2=seq(-6, 6, length.out = 100))
classes <- predict(svmfit_radial, test_df)
full_df <- test_df %>% 
  mutate(class=classes) %>% 
  mutate(type="simulated") %>% 
  bind_rows(x %>% mutate(type="actual"))
 
ggplot(data=full_df %>% filter(type=="simulated"),
       aes(x=V1, y=V2, colour=as.factor(class))) +
  geom_point(alpha=0.1) +
  geom_point(data=full_df %>% filter(type=="actual"))

Regression

  1. Suppose a linear regression model of the form:

\[\begin{equation} y_i \stackrel{i.i.d.}{\sim} \mathcal{N}(\alpha + \beta x_i, \sigma) \end{equation}\]

and assume you have access to \((y_i,x_i)_{i=1}^{K}\). Show that the maximum likelihood estimators of the parameters are the same as those minimising the mean squared loss function:

\[\begin{equation} L = \frac{1}{K} \sum_{i=1}^{K} (y_i - (\alpha + \beta x_i))^2 \end{equation}\]

The likelihood is:

\[\begin{equation} \mathcal{L} = \prod_{i=1}^{K} \frac{1}{\sqrt{2\pi\sigma^2}} \exp(-(y_i - (\alpha + \beta x_i))^2/2\sigma^2) \end{equation}\]

Writing down the log-likelihood, we have:

\[\begin{equation} \mathcal{l} = - \frac{K}{2}\log 2\pi - \frac{K}{2}\log \sigma^2 - \frac{1}{2\sigma^2} \sum_{i=1}^{K} (y_i - (\alpha + \beta x_i))^2 \end{equation}\]

We now note that the squared loss term in the log-likelihood is of the same form as that in the mean squared loss function, so the parameter estimators are the same in both cases.

  1. Generate simulated regression data via:

\[\begin{equation} y_i \stackrel{i.i.d.}{\sim} \mathcal{N}(\alpha + \beta x_i, \sigma) \end{equation}\]

where \(\alpha=1\), \(\beta=1\), \(\sigma=1\) and \(x_i \stackrel{i.i.d.}{\sim}\mathcal{N}(0, 4)\) for \(i =1,...,100\).

alpha <- 1
beta <- 1
sigma <- 1
n <- 100
x <- rnorm(n, 0, 4)
y <- alpha + beta * x + rnorm(n, 0, sigma)
plot(x, y)

  1. Create a function which implements one epoch of gradient descent and updates estimates of both \(\alpha\) and \(\beta\) _{i=1}^{K} (y_i - (+ x_i))
update_parameters <- function(alpha, beta, y, x, eta) {
  K <- length(y)
  dl_da <- -(2 / K) * sum(y - (alpha + beta * x))
  dl_db <- -(2 / K) * sum(x * (y - (alpha + beta * x)))
  
  alpha <- alpha - eta * dl_da
  beta <- beta - eta * dl_db
  return(c(alpha, beta))
}
  1. Update parameters across many epochs and visualise their parameter estimates during training
f_train_model <- function(n_epochs, eta, y, x, inits=c(0, 0)) {

  m_params <- matrix(nrow = n_epochs, ncol = 2)
  m_params[1, ] <- inits
  for(i in 2:n_epochs) {
    m_params[i, ] <- update_parameters(m_params[(i - 1), 1], m_params[(i - 1), 2], y, x, eta)
  }
  m_params <- as.data.frame(m_params)
  colnames(m_params) <- c("alpha", "beta")
  m_params <- m_params %>% 
    mutate(epoch=seq_along(alpha))
  return(m_params)
}
m_params <- f_train_model(2000, 0.01, y, x)
m_params %>% 
  melt(id.vars="epoch") %>% 
  ggplot(aes(x=epoch, y=value, colour=as.factor(variable))) +
  geom_line()

  1. Draw the loss surface and use it to visualise the path of parameter estimates during training
msl <- function(y, x, alpha, beta) {
  K <- length(y)
  return(1/K * sum((y-(alpha + beta * x))^2))
}

params <- expand_grid(alpha=seq(-2, 5, 0.1), beta=seq(-2, 5, 0.1))
z <- map2_dbl(params$alpha, params$beta, ~msl(y, x, .x, .y))
params <- params %>% 
  mutate(loss=z)

m_params <- m_params[1:100, ] %>% 
  mutate(type="training")

df1 <- params %>% 
  mutate(type="loss") %>% 
  bind_rows(m_params)

ggplot(df1 %>% filter(type=="loss"), aes(x=alpha, y=beta)) +
  geom_contour(aes(z=log(loss)), bins = 10) +
  geom_point(data=df1 %>% filter(type=="training")) +
  geom_line(data=df1 %>% filter(type=="training"))

  1. How does choice of \(\eta\) affect the rate of training?
# Low eta
f_train_model(2000, 0.00001, y, x) %>% 
  melt(id.vars="epoch") %>% 
  ggplot(aes(x=epoch, y=value, colour=as.factor(variable))) +
  geom_line()


# High eta
f_train_model(2000, 0.05, y, x) %>% 
  melt(id.vars="epoch") %>% 
  ggplot(aes(x=epoch, y=value, colour=as.factor(variable))) +
  geom_line()

The higher \(\eta\) is, the faster the rate of training.

  1. Why does a high value of \(\eta\) cause issues for training?
m_params <- f_train_model(100, 0.1, y, x)

m_params <- m_params[1:4, ] %>% 
  mutate(type="training")

df1 <- params %>% 
  mutate(type="loss") %>% 
  bind_rows(m_params)

ggplot(df1 %>% filter(type=="loss"), aes(x=alpha, y=beta)) +
  geom_contour(aes(z=log(loss)), bins = 10) +
  geom_point(data=df1 %>% filter(type=="training")) +
  geom_line(data=df1 %>% filter(type=="training"))

If too high a value is chosen, the steps taken in a given direction overshoot the curvature, resulting in an increase in the loss. If these steps are repeatedly taken, the loss divergences.

  1. For a reasonable value of \(\eta\), visualise the mean squared loss over time
m_params <- f_train_model(500, 0.01, y, x)

losses <- map_dbl(seq_along(m_params$alpha), ~msl(y, x, m_params$alpha[.], m_params$beta[.]))
m_params <- m_params %>% 
  mutate(loss=losses)
m_params %>% 
  ggplot(aes(x=epoch, y=loss)) +
  geom_line()

  1. Overlay your estimated regression line on the data
x_sim <- seq(-12, 12, 0.1)
y_sim <- m_params$alpha[nrow(m_params)] + m_params$beta[nrow(m_params)] * x_sim
df <- tibble(x, y, type="actual") %>% 
  bind_rows(tibble(x=x_sim, y=y_sim, type="simulated"))
ggplot(df %>% filter(type=="actual"), aes(x=x, y=y)) +
  geom_point() +
  geom_line(data=df %>% filter(type=="simulated"))

  1. Generate similar data except with \(\beta=0.02\) and \(x\sim\mathcal{N}(0, 400)\). How is training impacted for this model?
alpha <- 1
beta <- 0.02
sigma <- 10
n <- 100
x <- rnorm(n, 0, 400)
y <- alpha + beta * x + rnorm(n, 0, sigma)
plot(x, y)


params <- expand_grid(alpha=seq(-2, 5, 0.1), beta=seq(-0.05, 0.05, 0.005))
z <- map2_dbl(params$alpha, params$beta, ~msl(y, x, .x, .y))
params <- params %>% 
  mutate(loss=z)

m_params <- f_train_model(1000, 0.0000001, y, x) %>% 
  mutate(type="training")

df1 <- params %>% 
  mutate(type="loss") %>% 
  bind_rows(m_params)

ggplot(df1 %>% filter(type=="loss"), aes(x=alpha, y=beta)) +
  geom_contour(aes(z=log(loss)), bins = 10) +
  geom_point(data=df1 %>% filter(type=="training")) +
  geom_line(data=df1 %>% filter(type=="training"))

Training is very slow for \(\alpha\) due to the different parameter scales.

  1. Standardise both the \(y\) and \(x\) variables and return training. How is the rate of convergence now?
y_tilde <- scale(y)[, 1]
x_tilde <- scale(x)[, 1]

params <- expand_grid(alpha=seq(-2, 5, 0.1), beta=seq(-1, 2, 0.01))
z <- map2_dbl(params$alpha, params$beta, ~msl(y_tilde, x_tilde, .x, .y))
params <- params %>% 
  mutate(loss=z)

m_params <- f_train_model(1000, 0.01, y_tilde, x_tilde, inits = c(-1, -1)) %>% 
  mutate(type="training")

df1 <- params %>% 
  mutate(type="loss") %>% 
  bind_rows(m_params)

ggplot(df1 %>% filter(type=="loss"), aes(x=alpha, y=beta)) +
  geom_contour(aes(z=log(loss)), bins = 10) +
  geom_point(data=df1 %>% filter(type=="training")) +
  geom_line(data=df1 %>% filter(type=="training"))

Training now much faster: in summary, gradient descent much more efficient if data are on similar scales. So, feature scaling really important.

  1. Generate regression data via:

\[\begin{equation} y_i \stackrel{i.i.d.}{\sim} \mathcal{N}(\alpha + \beta x_i, \sigma) \end{equation}\]

where \(\alpha=1\), \(\beta=1\), \(\sigma=5\) and \(x_i \stackrel{i.i.d.}{\sim}\mathcal{N}(0, 4)\) for \(i =1,...,30\).

Use gradient descent to obtain parameter estimates for a quadratic regression and graph the fitted line versus the data.

alpha <- 1
beta <- 1
sigma <- 5
n <- 30
x <- rnorm(n, 0, 4)
y <- alpha + beta * x + rnorm(n, 0, sigma)

update_parameters_quad <- function(alpha, beta, gamma, y, x, eta) {
  K <- length(y)
  dl_da <- -(2 / K) * sum(y - (alpha + beta * x + gamma * x^2))
  dl_db <- -(2 / K) * sum(x * (y - (alpha + beta * x + gamma * x^2)))
  dl_dg <- -(2 / K) * sum(x^2 * (y - (alpha + beta * x + gamma * x^2)))
  
  alpha <- alpha - eta * dl_da
  beta <- beta - eta * dl_db
  gamma <- gamma - eta * dl_dg
  return(c(alpha, beta, gamma))
}

n_epochs <- 5000
m_params <- matrix(nrow = n_epochs, ncol = 3)
m_params[1, ] <- c(0, 0, 0)
for(i in 2:n_epochs) {
  m_params[i, ] <- update_parameters_quad(m_params[(i - 1), 1], m_params[(i - 1), 2],
                                          m_params[(i - 1), 3], y, x, 0.001)
}
m_params <- as.data.frame(m_params)
colnames(m_params) <- c("alpha", "beta", "gamma")
m_params <- m_params %>% 
  mutate(epoch=seq_along(alpha))
m_params %>% 
  melt(id.vars="epoch") %>% 
  ggplot(aes(x=epoch, y=value, colour=as.factor(variable))) +
  geom_line()


x_sim <- seq(-12, 12, 0.1)
y_sim <- m_params$alpha[nrow(m_params)] + m_params$beta[nrow(m_params)] * x_sim + m_params$gamma[nrow(m_params)] * x_sim^2

f_quadratic_regression <- function(x, alpha, beta, gamma) {
  return(alpha + beta * x + gamma * x^2)
}
a <- m_params$alpha[nrow(m_params)]
b <- m_params$beta[nrow(m_params)]
c <- m_params$gamma[nrow(m_params)]
f_quad <- function(x) f_quadratic_regression(x, a, b, c)

df <- tibble(x, y, type="actual") %>% 
  bind_rows(tibble(x=x_sim, y=y_sim, type="simulated"))
ggplot(df %>% filter(type=="actual"), aes(x=x, y=y)) +
  geom_point() +
  geom_line(data=df %>% filter(type=="simulated"))

  1. Fit the same data as generated in the previous question using linear regression. Compare the fit of the quadratic and linear models on new data generated from the same data generating process except:

\(x_i \stackrel{i.i.d.}{\sim}\mathcal{N}(5, 4)\) for \(i =1,...,30\)

Which has a better root mean squared accuracy on the new data? Why?

m_params <- f_train_model(2000, 0.01, y, x)

f_linear_regression <- function(x, alpha, beta) {
  return(alpha + beta * x)
}

f_linear <- function(x) f_linear_regression(x,m_params$alpha[nrow(m_params)],
                                             m_params$beta[nrow(m_params)])

x1 <- rnorm(n, 5, 4)
y1 <- alpha + beta * x1 + rnorm(n, 0, sigma)

# predict y and calculate RMSEs
y1_linear <- map_dbl(x1, f_linear)
y1_quad <- map_dbl(x1, f_quad)

sqrt(mean(y1_linear - y1)^2)
[1] 0.4445086
sqrt(mean(y1_quad - y1)^2)
[1] 3.150959

The linear model generalises better. The quadratic model has overfitted the training data, so doesn’t extrapolate well to new data regimes.

KNN regression

  1. Generate \(n=200\) paired \((x_i, y_i)\) data points by sampling:

\[\begin{equation} x_i \sim \mathcal{N}(0, 4) \end{equation}\]

and then:

\[\begin{equation} y_i \sim \mathcal{N}(sin(x_i), 0.2) \end{equation}\]

library(RANN)
library(tidyverse)
n <- 200
x <- rnorm(n, 0, 4)
y <- sin(x) + rnorm(n, 0, 0.2)

df <- tibble(x, y) %>% 
  mutate(type="actual")

ggplot(df %>% filter(type=="actual"), aes(x=x, y=y)) +
  geom_point()

  1. Create a function which returns the KNNs for a given \(x\) value (and \(k\))
# use RANN library (which implements KD trees)
x_tilde <- -2.5
fit <- nn2(x, c(x_tilde), k = 10)
idxs <- fit$nn.idx[1, ]

# visualise nearest neighbours around a point
tibble(x, y) %>% 
  mutate(idx=seq_along(x)) %>% 
  mutate(nn=idx%in%idxs) %>% 
  ggplot(aes(x=x, y=y, colour=as.factor(nn))) +
  geom_point()

  1. Using your answer to the previous code, create a KNN regression estimator using Euclidean distance with various values of \(k\) to predict \(y\). How does \(k\) influence the results?
f_kk_regression <- function(x_tilde, x, y, k=10) {
  fit <- nn2(x, c(x_tilde), k = k)
  idxs <- fit$nn.idx[1, ]
  return(mean(y[idxs]))
}

x_sim <- seq(-10, 10, 0.1)

# k = 1
y_sim <- map_dbl(x_sim, ~f_kk_regression(., x, y, k=10))

df <- tibble(x, y) %>% 
  mutate(type="actual") %>% 
  bind_rows(tibble(x=x_sim, y=y_sim) %>% mutate(type="regression"))

ggplot(df %>% filter(type=="actual"), aes(x=x, y=y)) +
  geom_point() +
  geom_line(data=df %>% filter(type=="regression"), colour="blue")


# k = 1
y_sim <- map_dbl(x_sim, ~f_kk_regression(., x, y, k=1))

df <- tibble(x, y) %>% 
  mutate(type="actual") %>% 
  bind_rows(tibble(x=x_sim, y=y_sim) %>% mutate(type="regression"))

ggplot(df %>% filter(type=="actual"), aes(x=x, y=y)) +
  geom_point() +
  geom_line(data=df %>% filter(type=="regression"), colour="blue")


# k = 100
y_sim <- map_dbl(x_sim, ~f_kk_regression(., x, y, k=100))

df <- tibble(x, y) %>% 
  mutate(type="actual") %>% 
  bind_rows(tibble(x=x_sim, y=y_sim) %>% mutate(type="regression"))

ggplot(df %>% filter(type=="actual"), aes(x=x, y=y)) +
  geom_point() +
  geom_line(data=df %>% filter(type=="regression"), colour="blue")

LS0tCnRpdGxlOiAiU3VwZXJ2aXNlZCBNTCBwcm9ibGVtcyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CnJtKGxpc3Q9bHMoKSkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkoZTEwNzEpCmBgYAoKCiMgU3VwcG9ydCB2ZWN0b3IgbWFjaGluZXMKMS4gU3VwcG9zZSBhIGNhdGVnb3JpY2FsIGxhYmVsOiAkeVxpblx7MCwxXH0kIGFuZCBhIGRhdGEgZ2VuZXJhdGluZyBwcm9jZXNzOiAkeHx5PTBcc2ltXG1hdGhjYWx7Tn0oMCwxKSQgYW5kICR4fHk9MVxzaW1cbWF0aGNhbHtOfSgyLCAxKSQuIEZpbmQgdGhlIG9wdGltYWwgY2xhc3NpZmljYXRpb24gYm91bmRhcnkuIEluIGRvaW5nIHNvLCBkZXRlcm1pbmUgdGhlIEJheWVzIG9wdGltYWwgY2xhc3NpZmljYXRpb24gYWNjdXJhY3kuIEluIGJvdGggY2FzZXMsIGFzc3VtZSBjYXRlZ29yaWVzIGFyZSBlcXVhbGx5IGFzIGxpa2VseSB0byBvY2N1ciBpbiB0aGUgZGF0YS4KCkRyYXdpbmcgdGhlIHR3byBkZW5zaXRpZXMKYGBge3J9Cm4gPC0gMTAwCnggPC0gc2VxKC0yLCA1LCBsZW5ndGgub3V0ID0gbikKeTAgPC0gbWFwX2RibCh4LCB+ZG5vcm0oLikpCnkxIDwtIG1hcF9kYmwoeCwgfmRub3JtKC4sIDIsIDEpKQp0aWJibGUoeCwgeTAsIHkxKSAlPiUgCiAgbWVsdChpZC52YXJzPSJ4IikgJT4lIAogIGdncGxvdChhZXMoeCwgeT12YWx1ZSwgY29sb3VyPWFzLmZhY3Rvcih2YXJpYWJsZSkpKSArCiAgZ2VvbV9saW5lKCkKYGBgCgpXZSBzZWUgdGhhdCB0aGUgdGhlIGRlbnNpdGllcyBjcm9zc292ZXIgYXQgYSBwb2ludC4gTGV0J3MgZmluZCB3aGF0IGl0IGlzLgoKJFxtYXRoY2Fse059KHh8MCwxKSA9IFxtYXRoY2Fse059KHh8MiwxKSQKCndoaWNoLCBieSBzeW1tZXRyeSwgb2NjdXJzIHdoZW4gJHg9MSQuCgpTbyB0aGUgb3B0aW1hbCBib3VuZGFyeSBpcyBwcmVkaWN0ICR5PTAkIHdoZW4gJHggPCAxJCBhbmQgJHk9MSQgd2hlbiAkeD4xJC4gQXQgdGhlIGJvdW5kYXJ5LCBwcmVkaWN0IGVpdGhlciB3aXRoIHByb2JhYmlsaXR5ICQxLzIkLgoKVGhlIG9wdGltYWwgY2xhc3NpZmljYXRpb24gYWNjdXJhY3kgaXMgJFByKHg8MXx5PTApID0gXGludF97LVxpbmZ0eX1eezF9IFxtYXRoY2Fse059KHh8MCwxKSBcYXBwcm94IDg0XCUkLgoKMi4gQXNzdW1lIHRoYXQgJCh4XzEseF8yKSd8eT0wXHNpbVxtYXRoY2Fse059KCgwLDApJywgSSkkIGFuZCAkKHhfMSx4XzIpJ3x5PTFcc2ltXG1hdGhjYWx7Tn0oKDEsIDEpJyxJKSQuIFNpbXVsYXRlIDEwMCBkcmF3cyBmcm9tIGVhY2ggY2F0ZWdvcnkuIFVzZSB0aGUgcGFyYW1ldHJpYyBmb3JtIG9mIGVhY2ggY29uZGl0aW9uYWwgZGVuc2l0eSB0byBjbGFzc2lmeSBlYWNoIHBvaW50IGFuZCB2aXN1YWxpc2UgdGhlc2U7IGluIGRvaW5nIHNvLCBjYWxjdWxhdGUgdGhlIG9wdGltYWwgY2xhc3NpZmljYXRpb24gcmF0ZSBvbiB0aGVzZSBkYXRhLgoKYGBge3J9CmxpYnJhcnkobXZ0bm9ybSkKcm12cm5vcm0yRCA8LSBmdW5jdGlvbihuLCBtdXgsIG11eSwgc2lnbWF4LCBzaWdtYXksIHJobyl7CiAgcmV0dXJuKHJtdm5vcm0obiwgYyhtdXgsIG11eSksCiAgICAgICAgICAgICAgICAgbWF0cml4KGMoc2lnbWF4XjIsIHNpZ21heCAqIHNpZ21heSAqIHJobywKICAgICAgICAgICAgICAgICAgICAgICAgICBzaWdtYXggKiBzaWdtYXkgKiByaG8sIHNpZ21heV4yKSwKICAgICAgICAgICAgICAgICAgICAgICAgbmNvbCA9IDIpKSkKfQoKbiA8LSAxMDAKeDAgPC0gcm12cm5vcm0yRChuLCAwLCAwLCAxLCAxLCAwKSAlPiUgCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICBtdXRhdGUoY2xhc3M9MCkKeDEgPC0gcm12cm5vcm0yRChuLCAxLCAxLCAxLCAxLCAwKSAlPiUgCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICBtdXRhdGUoY2xhc3M9MSkKCiMgdHJ1ZSBjbGFzc2VzCnggPC0geDAgJT4lIAogIGJpbmRfcm93cyh4MSkKeCAlPiUgCiAgZ2dwbG90KGFlcyh4PVYxLCB5PVYyLCBjb2xvdXI9YXMuZmFjdG9yKGNsYXNzKSkpICsKICBnZW9tX3BvaW50KCkKYGBgCgpDYWxjdWxhdGUgdGhlIHByb2JhYmlsaXR5IGRlbnNpdHkgZm9yIGVhY2ggcG9pbnQuIEluIGRvaW5nIHNvLCB0aGlzIHlpZWxkcyBhIHByb2JhYmlsaXR5IG9mIGJlbG9uZ2luZyB0byBlYWNoIGNsYXNzICh3aGVyZSB3ZSBrbm93IGVhY2ggYXJlIGVxdWFsbHkgbGlrZWx5KToKCiRQcih5PTF8eF8xLHhfMikgPSBcZnJhY3twKHhfMSx4XzJ8eT0xKSBQcih5PTEpfXtwKHhfMSx4XzIpfT1cZnJhY3twKHhfMSx4XzJ8eT0xKSBQcih5PTEpfXtwKHhfMSx4XzJ8eT0xKSBQcih5PTEpICsgcCh4XzEseF8yfHk9MCkgUHIoeT0wKX0kCgpgYGB7cn0KZG12cm5vcm0yRCA8LSBmdW5jdGlvbih4LCBtdXgsIG11eSwgc2lnbWF4LCBzaWdtYXksIHJobywgbG9nPUZBTFNFKXsKICByZXR1cm4oZG12bm9ybSh4LCBjKG11eCwgbXV5KSwKICAgICAgICAgICAgICAgICBtYXRyaXgoYyhzaWdtYXheMiwgc2lnbWF4ICogc2lnbWF5ICogcmhvLAogICAgICAgICAgICAgICAgICAgICAgICAgIHNpZ21heCAqIHNpZ21heSAqIHJobywgc2lnbWF5XjIpLAogICAgICAgICAgICAgICAgICAgICAgICAgbmNvbCA9IDIpLAogICAgICAgICAgICAgICAgIGxvZykpCn0KCnAwIDwtIG1hcF9kYmwoc2VxX2Fsb25nKHgkVjEpLCB+ZG12cm5vcm0yRCh4Wy4sIDE6Ml0sIDAsIDAsIDEsIDEsIDApKQpwMSA8LSBtYXBfZGJsKHNlcV9hbG9uZyh4JFYxKSwgfmRtdnJub3JtMkQoeFsuLCAxOjJdLCAxLCAxLCAxLCAxLCAwKSkKCmRmIDwtIHggJT4lIAogIG11dGF0ZShwMD1wMCwgcDE9cDEpICU+JSAKICBtdXRhdGUoY2xhc3NfcHJlZD1pZl9lbHNlKHAxID4gcDAsIDEsIDApKSAlPiUgCiAgbXV0YXRlKGNvcnJlY3Q9aWZfZWxzZShjbGFzcz09Y2xhc3NfcHJlZCwgMSwgMCkpCgptZWFuKGRmJGNvcnJlY3QpCmBgYAoKVmlzdWFsaXNlIGluY29ycmVjdCBwcmVkaWN0aW9ucwpgYGB7cn0KZGYgJT4lIAogIGdncGxvdChhZXMoeD1WMSwgeT1WMiwgY29sb3VyPWFzLmZhY3Rvcihjb3JyZWN0KSkpICsKICBnZW9tX3BvaW50KCkKYGBgCgozLiBXaGF0IGlzIHRoZSBvcHRpbWFsIGNsYXNzaWZpY2F0aW9uIGJvdW5kYXJ5IGhlcmU/CgpUaGUgb3B0aW1hbCBjbGFzc2lmaWNhdGlvbiBib3VuZGFyeSBpcyBnaXZlbiBieToKClxiZWdpbntlcXVhdGlvbn0KXG1hdGhjYWx7Tn0oeF8xLHhfMnwoMCwwKScsIEkpID0gXG1hdGhjYWx7Tn0oeF8xLHhfMnwoMSwxKScsIEkpClxlbmR7ZXF1YXRpb259CgpTb2x2aW5nIHRoZXNlIGVxdWF0aW9ucywgd2UgZmluZCB0aGF0OiAkeF8yPSAxIC0geF8xJC4gUGxvdHRpbmcgdGhpcyBvbiB0aGUgZGF0YSBkaXN0cmlidXRpb246CmBgYHtyfQp4ICU+JSAKICBnZ3Bsb3QoYWVzKHg9VjEsIHk9VjIsIGNvbG91cj1hcy5mYWN0b3IoY2xhc3MpKSkgKwogIGdlb21fcG9pbnQoKSArCiAgc3RhdF9mdW5jdGlvbihmdW4gPSBmdW5jdGlvbihWMSkgMSAtIFYxLCBjb2xvdXI9ImJsYWNrIikgKwogIGNvb3JkX2ZpeGVkKCkgKwogIHNjYWxlX3lfY29udGludW91cyhsaW1pdHMgPSBjKC00LCA0KSkgKwogIHNjYWxlX3hfY29udGludW91cyhsaW1pdHMgPSBjKC00LCA0KSkKYGBgCgo0LiBBcHBseSAobGluZWFyKSBTVk0gdG8gdGhlIGRhdGFzZXQgeW91IGdlbmVyYXRlZCBpbiB0aGUgcHJldmlvdXMgcXVlc3Rpb24uIFdoYXQncyBpdHMgcHJlZGljdGl2ZSBhY2N1cmFjeSBvbiB5b3VyIGRhdGFzZXQ/IEhvdyBjbG9zZSBkb2VzIGl0IGdldCB0byB0aGUgb3B0aW1hbCBjbGFzc2lmaWNhdGlvbiBib3VuZGFyeT8KCk5vdCBiYWQgYnV0IG92ZXJmaXRzIGEgYml0LgpgYGB7cn0KZGYgPC0gZGYgJT4lIAogIG11dGF0ZShjbGFzcz1hcy5mYWN0b3IoY2xhc3MpKQpzdm1maXQgPC0gc3ZtKGNsYXNzIH4gLiwgZGF0YSA9IGRmICU+JSBzZWxlY3QoVjEsIFYyLCBjbGFzcyksIGtlcm5lbCA9ICJsaW5lYXIiLCBjb3N0ID0gMTAsIHNjYWxlID0gRkFMU0UsCiAgICAgICAgICAgICAgcHJvYmFiaWxpdHk9VFJVRSkKCiMjIHByZWRpY3RpdmUgYWNjdXJhY3kKY2xhc3MxIDwtIHByZWRpY3Qoc3ZtZml0LCBkZikKZGYgJT4lIAogIHNlbGVjdChjbGFzcykgJT4lIAogIG11dGF0ZShjbGFzc19wcmVkPWNsYXNzMSkgJT4lIAogIG11dGF0ZShjb3JyZWN0PWlmX2Vsc2UoY2xhc3M9PWNsYXNzX3ByZWQsIDEsIDApKSAlPiUgCiAgc3VtbWFyaXNlKG1lYW4oY29ycmVjdCkpCgojIHBsb3QgY2xhc3NpZmljYXRpb24gYm91bmRhcnkKdGVzdF9kZiA8LSBleHBhbmRfZ3JpZChWMT1zZXEoLTQsIDQsIGxlbmd0aC5vdXQgPSAxMDApLCBWMj1zZXEoLTQsIDQsIGxlbmd0aC5vdXQgPSAxMDApKQpjbGFzc2VzIDwtIHByZWRpY3Qoc3ZtZml0LCB0ZXN0X2RmKQp0ZXN0X2RmICU+JSAKICBtdXRhdGUoY2xhc3NfcHJlZD1jbGFzc2VzKSAlPiUgCiAgZ2dwbG90KGFlcyh4PVYxLCB5PVYyLCBjb2xvdXI9YXMuZmFjdG9yKGNsYXNzX3ByZWQpKSkgKwogIGdlb21fcG9pbnQoKSArCiAgc3RhdF9mdW5jdGlvbihmdW4gPSBmdW5jdGlvbihWMSkgMSAtIFYxLCBjb2xvdXI9ImJsYWNrIikgKwogIGNvb3JkX2ZpeGVkKCkgKwogIHNjYWxlX3lfY29udGludW91cyhsaW1pdHMgPSBjKC00LCA0KSkgKwogIHNjYWxlX3hfY29udGludW91cyhsaW1pdHMgPSBjKC00LCA0KSkKYGBgCgpIb3cgZG8gcHJlZGljdGVkIGNsYXNzIHByb2JhYmlsaXRpZXMgY29tcGFyZSB3aXRoIHRydWUgb25lcz8KYGBge3J9CmNsYXNzX3Byb2JzIDwtIHByZWRpY3Qoc3ZtZml0LCBkZiwgcHJvYmFiaWxpdHkgPSBUKQpgYGAKCgoKNS4gRml0IGFuIFNWTSB3aXRoIGEgcmFkaWFsIGJhc2lzIGZ1bmN0aW9uIGtlcm5lbC4gV2hhdCdzIGl0cyBwcmVkaWN0aXZlIGFjY3VyYWN5PyBXaGF0IGRvIGl0cyBjbGFzc2lmaWNhdGlvbiBib3VuZGFyaWVzIGxvb2sgbGlrZT8KYGBge3J9CmRmIDwtIGRmICU+JSAKICBtdXRhdGUoY2xhc3M9YXMuZmFjdG9yKGNsYXNzKSkKc3ZtZml0X3JhZGlhbCA8LSBzdm0oY2xhc3MgfiAuLCBkYXRhID0gZGYgJT4lIHNlbGVjdChWMSwgVjIsIGNsYXNzKSwga2VybmVsID0gInJhZGlhbCIsIGNvc3QgPSAxMCwgc2NhbGUgPSBGQUxTRSkKCiMjIHByZWRpY3RpdmUgYWNjdXJhY3kKY2xhc3MxIDwtIHByZWRpY3Qoc3ZtZml0X3JhZGlhbCwgZGYpCmRmICU+JSAKICBzZWxlY3QoY2xhc3MpICU+JSAKICBtdXRhdGUoY2xhc3NfcHJlZD1jbGFzczEpICU+JSAKICBtdXRhdGUoY29ycmVjdD1pZl9lbHNlKGNsYXNzPT1jbGFzc19wcmVkLCAxLCAwKSkgJT4lIAogIHN1bW1hcmlzZShtZWFuKGNvcnJlY3QpKQoKIyBwbG90IGNsYXNzaWZpY2F0aW9uIGJvdW5kYXJ5CnRlc3RfZGYgPC0gZXhwYW5kX2dyaWQoVjE9c2VxKC00LCA0LCBsZW5ndGgub3V0ID0gMTAwKSwgVjI9c2VxKC00LCA0LCBsZW5ndGgub3V0ID0gMTAwKSkKY2xhc3NlcyA8LSBwcmVkaWN0KHN2bWZpdF9yYWRpYWwsIHRlc3RfZGYpCnRlc3RfZGYgJT4lIAogIG11dGF0ZShjbGFzc19wcmVkPWNsYXNzZXMpICU+JSAKICBnZ3Bsb3QoYWVzKHg9VjEsIHk9VjIsIGNvbG91cj1hcy5mYWN0b3IoY2xhc3NfcHJlZCkpKSArCiAgZ2VvbV9wb2ludCgpICsKICBzdGF0X2Z1bmN0aW9uKGZ1biA9IGZ1bmN0aW9uKFYxKSAxIC0gVjEsIGNvbG91cj0iYmxhY2siKSArCiAgY29vcmRfZml4ZWQoKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGxpbWl0cyA9IGMoLTQsIDQpKSArCiAgc2NhbGVfeF9jb250aW51b3VzKGxpbWl0cyA9IGMoLTQsIDQpKQpgYGAKCjYuIENvbXBhcmUgdGhlIHByZWRpY3RpdmUgYWNjdXJhY3kgb2YgZWFjaCBTVk0gbW9kZWwgb24gYSBuZXcgZGF0YXNldCBkcmF3biBmcm9tIHRoZSBzYW1lIGRpc3RyaWJ1dGlvbi4gV2hpY2ggU1ZNIHBlcmZvcm1zIGJlc3Q/CgpJdCBzaG91bGQgYmUgdGhlIGxpbmVhciBvbmUgYXMgdGhlIHJhZGlhbCBvbmUgb3ZlcmZpdHMgdGhlIGRhdGEuCmBgYHtyfQpuIDwtIDEwMAp4MCA8LSBybXZybm9ybTJEKG4sIDAsIDAsIDEsIDEsIDApICU+JSAKICBhcy5kYXRhLmZyYW1lKCkgJT4lIAogIG11dGF0ZShjbGFzcz0wKQp4MSA8LSBybXZybm9ybTJEKG4sIDEsIDEsIDEsIDEsIDApICU+JSAKICBhcy5kYXRhLmZyYW1lKCkgJT4lIAogIG11dGF0ZShjbGFzcz0xKQp4IDwtIHgwICU+JSAKICBiaW5kX3Jvd3MoeDEpCmNsYXNzX2xpbmVhciA8LSBwcmVkaWN0KHN2bWZpdCwgeCkKY2xhc3NfcmFkaWFsIDwtIHByZWRpY3Qoc3ZtZml0X3JhZGlhbCwgeCkKCnggJT4lIAogIG11dGF0ZShjbGFzc19saW5lYXI9Y2xhc3NfbGluZWFyLAogICAgICAgICBjbGFzc19yYWRpYWw9Y2xhc3NfcmFkaWFsKSAlPiUgCiAgc3VtbWFyaXNlKG1lYW4oY2xhc3M9PWNsYXNzX2xpbmVhciksIG1lYW4oY2xhc3M9PWNsYXNzX3JhZGlhbCkpCmBgYAoKNy4gR2VuZXJhdGUgZGF0YSBmb3Igb25lIGNsYXNzIGJ5IHNhbXBsaW5nICQoeF8xLHhfMiknfHk9MFxzaW1cbWF0aGNhbHtOfSgoMCwwKScsIEkpJC4gRm9yIHRoZSBvdGhlciBjbGFzcyBzYW1wbGUgZnJvbSAkKHhfMSx4XzIpJ3x5PTFcc2ltXG1hdGhjYWx7Tn0oKDAsMCknLCAyIEkpJCBzdWNoIHRoYXQgJHhfMV4yICsgeF8yXjIgXGdlcSAzXjIkLgpgYGB7cn0KbiA8LSAxMDAKeDAgPC0gcm12cm5vcm0yRChuLCAwLCAwLCAxLCAxLCAwKSAlPiUgCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICBtdXRhdGUoY2xhc3M9MCkKCmZfc2FtcGxlX2NvbnN0cmFpbmVkIDwtIGZ1bmN0aW9uKHIpIHsKICB4X3Rlc3QgPC0gcm12cm5vcm0yRCgxLCAwLCAwLCAyLCAyLCAwKVsxLCBdCiAgd2hpbGUoKHhfdGVzdFsxXV4yICsgeF90ZXN0WzJdXjIpIDwgcl4yKQogICAgeF90ZXN0IDwtIHJtdnJub3JtMkQoMSwgMCwgMCwgMiwgMiwgMClbMSwgXQogIHJldHVybih4X3Rlc3QpCn0KCngxIDwtIG1hdHJpeChucm93ID0gbiwgbmNvbCA9IDIpCmZvcihpIGluIDE6bikKICB4MVtpLCBdIDwtIGZfc2FtcGxlX2NvbnN0cmFpbmVkKDMpCgp4MSA8LSB4MSAlPiUgCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICBtdXRhdGUoY2xhc3M9MSkKCnggPC0geDAgJT4lIAogIGJpbmRfcm93cyh4MSkgJT4lIAogIG11dGF0ZShjbGFzcz1hcy5mYWN0b3IoY2xhc3MpKQp4ICU+JSAKICBnZ3Bsb3QoYWVzKHg9VjEsIHk9VjIsIGNvbG91cj1jbGFzcykpICsKICBnZW9tX3BvaW50KCkKYGBgCgo4LiBGaXQgYSBjbGFzc2lmaWVyIHVzaW5nIHJhZGlhbCBiYXNpcyBmdW5jdGlvbnMgYW5kIHZpc3VhbGlzZSB0aGUgZGVjaXNpb24gYm91bmRhcnkKYGBge3J9CnN2bWZpdF9yYWRpYWwgPC0gc3ZtKGNsYXNzIH4gLiwgZGF0YSA9IHgsIGtlcm5lbCA9ICJyYWRpYWwiLCBjb3N0ID0gMTAsIHNjYWxlID0gRkFMU0UpCgp0ZXN0X2RmIDwtIGV4cGFuZF9ncmlkKFYxPXNlcSgtNiwgNiwgbGVuZ3RoLm91dCA9IDEwMCksIFYyPXNlcSgtNiwgNiwgbGVuZ3RoLm91dCA9IDEwMCkpCmNsYXNzZXMgPC0gcHJlZGljdChzdm1maXRfcmFkaWFsLCB0ZXN0X2RmKQpmdWxsX2RmIDwtIHRlc3RfZGYgJT4lIAogIG11dGF0ZShjbGFzcz1jbGFzc2VzKSAlPiUgCiAgbXV0YXRlKHR5cGU9InNpbXVsYXRlZCIpICU+JSAKICBiaW5kX3Jvd3MoeCAlPiUgbXV0YXRlKHR5cGU9ImFjdHVhbCIpKQogCmdncGxvdChkYXRhPWZ1bGxfZGYgJT4lIGZpbHRlcih0eXBlPT0ic2ltdWxhdGVkIiksCiAgICAgICBhZXMoeD1WMSwgeT1WMiwgY29sb3VyPWFzLmZhY3RvcihjbGFzcykpKSArCiAgZ2VvbV9wb2ludChhbHBoYT0wLjEpICsKICBnZW9tX3BvaW50KGRhdGE9ZnVsbF9kZiAlPiUgZmlsdGVyKHR5cGU9PSJhY3R1YWwiKSkKYGBgCgojIFJlZ3Jlc3Npb24KCjEuIFN1cHBvc2UgYSBsaW5lYXIgcmVncmVzc2lvbiBtb2RlbCBvZiB0aGUgZm9ybToKClxiZWdpbntlcXVhdGlvbn0KeV9pIFxzdGFja3JlbHtpLmkuZC59e1xzaW19IFxtYXRoY2Fse059KFxhbHBoYSArIFxiZXRhIHhfaSwgXHNpZ21hKQpcZW5ke2VxdWF0aW9ufQoKYW5kIGFzc3VtZSB5b3UgaGF2ZSBhY2Nlc3MgdG8gJCh5X2kseF9pKV97aT0xfV57S30kLiBTaG93IHRoYXQgdGhlIG1heGltdW0gbGlrZWxpaG9vZCBlc3RpbWF0b3JzIG9mIHRoZSBwYXJhbWV0ZXJzIGFyZSB0aGUgc2FtZSBhcyB0aG9zZSBtaW5pbWlzaW5nIHRoZSBtZWFuIHNxdWFyZWQgbG9zcyBmdW5jdGlvbjoKClxiZWdpbntlcXVhdGlvbn0KTCA9IFxmcmFjezF9e0t9IFxzdW1fe2k9MX1ee0t9ICh5X2kgLSAoXGFscGhhICsgXGJldGEgeF9pKSleMgpcZW5ke2VxdWF0aW9ufQoKVGhlIGxpa2VsaWhvb2QgaXM6CgpcYmVnaW57ZXF1YXRpb259ClxtYXRoY2Fse0x9ID0gXHByb2Rfe2k9MX1ee0t9IFxmcmFjezF9e1xzcXJ0ezJccGlcc2lnbWFeMn19IFxleHAoLSh5X2kgLSAoXGFscGhhICsgXGJldGEgeF9pKSleMi8yXHNpZ21hXjIpClxlbmR7ZXF1YXRpb259CgpXcml0aW5nIGRvd24gdGhlIGxvZy1saWtlbGlob29kLCB3ZSBoYXZlOgoKXGJlZ2lue2VxdWF0aW9ufQpcbWF0aGNhbHtsfSA9IC0gXGZyYWN7S317Mn1cbG9nIDJccGkgLSBcZnJhY3tLfXsyfVxsb2cgXHNpZ21hXjIgLSBcZnJhY3sxfXsyXHNpZ21hXjJ9IFxzdW1fe2k9MX1ee0t9ICh5X2kgLSAoXGFscGhhICsgXGJldGEgeF9pKSleMgpcZW5ke2VxdWF0aW9ufQoKV2Ugbm93IG5vdGUgdGhhdCB0aGUgc3F1YXJlZCBsb3NzIHRlcm0gaW4gdGhlIGxvZy1saWtlbGlob29kIGlzIG9mIHRoZSBzYW1lIGZvcm0gYXMgdGhhdCBpbiB0aGUgbWVhbiBzcXVhcmVkIGxvc3MgZnVuY3Rpb24sIHNvIHRoZSBwYXJhbWV0ZXIgZXN0aW1hdG9ycyBhcmUgdGhlIHNhbWUgaW4gYm90aCBjYXNlcy4KCgoyLiBHZW5lcmF0ZSBzaW11bGF0ZWQgcmVncmVzc2lvbiBkYXRhIHZpYToKClxiZWdpbntlcXVhdGlvbn0KeV9pIFxzdGFja3JlbHtpLmkuZC59e1xzaW19IFxtYXRoY2Fse059KFxhbHBoYSArIFxiZXRhIHhfaSwgXHNpZ21hKQpcZW5ke2VxdWF0aW9ufQoKd2hlcmUgJFxhbHBoYT0xJCwgJFxiZXRhPTEkLCAkXHNpZ21hPTEkIGFuZCAkeF9pIFxzdGFja3JlbHtpLmkuZC59e1xzaW19XG1hdGhjYWx7Tn0oMCwgNCkkIGZvciAkaSA9MSwuLi4sMTAwJC4KCmBgYHtyfQphbHBoYSA8LSAxCmJldGEgPC0gMQpzaWdtYSA8LSAxCm4gPC0gMTAwCnggPC0gcm5vcm0obiwgMCwgNCkKeSA8LSBhbHBoYSArIGJldGEgKiB4ICsgcm5vcm0obiwgMCwgc2lnbWEpCnBsb3QoeCwgeSkKYGBgCgozLiBDcmVhdGUgYSBmdW5jdGlvbiB3aGljaCBpbXBsZW1lbnRzIG9uZSBlcG9jaCBvZiBncmFkaWVudCBkZXNjZW50IGFuZCB1cGRhdGVzIGVzdGltYXRlcyBvZiBib3RoICRcYWxwaGEkIGFuZCAkXGJldGEkClxmcmFjezJ9e0t9XHN1bV97aT0xfV57S30gKHlfaSAtIChcYWxwaGEgKyBcYmV0YSB4X2kpKQpgYGB7cn0KdXBkYXRlX3BhcmFtZXRlcnMgPC0gZnVuY3Rpb24oYWxwaGEsIGJldGEsIHksIHgsIGV0YSkgewogIEsgPC0gbGVuZ3RoKHkpCiAgZGxfZGEgPC0gLSgyIC8gSykgKiBzdW0oeSAtIChhbHBoYSArIGJldGEgKiB4KSkKICBkbF9kYiA8LSAtKDIgLyBLKSAqIHN1bSh4ICogKHkgLSAoYWxwaGEgKyBiZXRhICogeCkpKQogIAogIGFscGhhIDwtIGFscGhhIC0gZXRhICogZGxfZGEKICBiZXRhIDwtIGJldGEgLSBldGEgKiBkbF9kYgogIHJldHVybihjKGFscGhhLCBiZXRhKSkKfQpgYGAKCjQuIFVwZGF0ZSBwYXJhbWV0ZXJzIGFjcm9zcyBtYW55IGVwb2NocyBhbmQgdmlzdWFsaXNlIHRoZWlyIHBhcmFtZXRlciBlc3RpbWF0ZXMgZHVyaW5nIHRyYWluaW5nCmBgYHtyfQpmX3RyYWluX21vZGVsIDwtIGZ1bmN0aW9uKG5fZXBvY2hzLCBldGEsIHksIHgsIGluaXRzPWMoMCwgMCkpIHsKCiAgbV9wYXJhbXMgPC0gbWF0cml4KG5yb3cgPSBuX2Vwb2NocywgbmNvbCA9IDIpCiAgbV9wYXJhbXNbMSwgXSA8LSBpbml0cwogIGZvcihpIGluIDI6bl9lcG9jaHMpIHsKICAgIG1fcGFyYW1zW2ksIF0gPC0gdXBkYXRlX3BhcmFtZXRlcnMobV9wYXJhbXNbKGkgLSAxKSwgMV0sIG1fcGFyYW1zWyhpIC0gMSksIDJdLCB5LCB4LCBldGEpCiAgfQogIG1fcGFyYW1zIDwtIGFzLmRhdGEuZnJhbWUobV9wYXJhbXMpCiAgY29sbmFtZXMobV9wYXJhbXMpIDwtIGMoImFscGhhIiwgImJldGEiKQogIG1fcGFyYW1zIDwtIG1fcGFyYW1zICU+JSAKICAgIG11dGF0ZShlcG9jaD1zZXFfYWxvbmcoYWxwaGEpKQogIHJldHVybihtX3BhcmFtcykKfQptX3BhcmFtcyA8LSBmX3RyYWluX21vZGVsKDIwMDAsIDAuMDEsIHksIHgpCm1fcGFyYW1zICU+JSAKICBtZWx0KGlkLnZhcnM9ImVwb2NoIikgJT4lIAogIGdncGxvdChhZXMoeD1lcG9jaCwgeT12YWx1ZSwgY29sb3VyPWFzLmZhY3Rvcih2YXJpYWJsZSkpKSArCiAgZ2VvbV9saW5lKCkKYGBgCgo1LiBEcmF3IHRoZSBsb3NzIHN1cmZhY2UgYW5kIHVzZSBpdCB0byB2aXN1YWxpc2UgdGhlIHBhdGggb2YgcGFyYW1ldGVyIGVzdGltYXRlcyBkdXJpbmcgdHJhaW5pbmcKYGBge3J9Cm1zbCA8LSBmdW5jdGlvbih5LCB4LCBhbHBoYSwgYmV0YSkgewogIEsgPC0gbGVuZ3RoKHkpCiAgcmV0dXJuKDEvSyAqIHN1bSgoeS0oYWxwaGEgKyBiZXRhICogeCkpXjIpKQp9CgpwYXJhbXMgPC0gZXhwYW5kX2dyaWQoYWxwaGE9c2VxKC0yLCA1LCAwLjEpLCBiZXRhPXNlcSgtMiwgNSwgMC4xKSkKeiA8LSBtYXAyX2RibChwYXJhbXMkYWxwaGEsIHBhcmFtcyRiZXRhLCB+bXNsKHksIHgsIC54LCAueSkpCnBhcmFtcyA8LSBwYXJhbXMgJT4lIAogIG11dGF0ZShsb3NzPXopCgptX3BhcmFtcyA8LSBtX3BhcmFtc1sxOjEwMCwgXSAlPiUgCiAgbXV0YXRlKHR5cGU9InRyYWluaW5nIikKCmRmMSA8LSBwYXJhbXMgJT4lIAogIG11dGF0ZSh0eXBlPSJsb3NzIikgJT4lIAogIGJpbmRfcm93cyhtX3BhcmFtcykKCmdncGxvdChkZjEgJT4lIGZpbHRlcih0eXBlPT0ibG9zcyIpLCBhZXMoeD1hbHBoYSwgeT1iZXRhKSkgKwogIGdlb21fY29udG91cihhZXMoej1sb2cobG9zcykpLCBiaW5zID0gMTApICsKICBnZW9tX3BvaW50KGRhdGE9ZGYxICU+JSBmaWx0ZXIodHlwZT09InRyYWluaW5nIikpICsKICBnZW9tX2xpbmUoZGF0YT1kZjEgJT4lIGZpbHRlcih0eXBlPT0idHJhaW5pbmciKSkKYGBgCgo2LiBIb3cgZG9lcyBjaG9pY2Ugb2YgJFxldGEkIGFmZmVjdCB0aGUgcmF0ZSBvZiB0cmFpbmluZz8KYGBge3J9CiMgTG93IGV0YQpmX3RyYWluX21vZGVsKDIwMDAsIDAuMDAwMDEsIHksIHgpICU+JSAKICBtZWx0KGlkLnZhcnM9ImVwb2NoIikgJT4lIAogIGdncGxvdChhZXMoeD1lcG9jaCwgeT12YWx1ZSwgY29sb3VyPWFzLmZhY3Rvcih2YXJpYWJsZSkpKSArCiAgZ2VvbV9saW5lKCkKCiMgSGlnaCBldGEKZl90cmFpbl9tb2RlbCgyMDAwLCAwLjA1LCB5LCB4KSAlPiUgCiAgbWVsdChpZC52YXJzPSJlcG9jaCIpICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZXBvY2gsIHk9dmFsdWUsIGNvbG91cj1hcy5mYWN0b3IodmFyaWFibGUpKSkgKwogIGdlb21fbGluZSgpCmBgYAoKVGhlIGhpZ2hlciAkXGV0YSQgaXMsIHRoZSBmYXN0ZXIgdGhlIHJhdGUgb2YgdHJhaW5pbmcuCgo3LiBXaHkgZG9lcyBhIGhpZ2ggdmFsdWUgb2YgJFxldGEkIGNhdXNlIGlzc3VlcyBmb3IgdHJhaW5pbmc/CgpgYGB7cn0KbV9wYXJhbXMgPC0gZl90cmFpbl9tb2RlbCgxMDAsIDAuMSwgeSwgeCkKCm1fcGFyYW1zIDwtIG1fcGFyYW1zWzE6NCwgXSAlPiUgCiAgbXV0YXRlKHR5cGU9InRyYWluaW5nIikKCmRmMSA8LSBwYXJhbXMgJT4lIAogIG11dGF0ZSh0eXBlPSJsb3NzIikgJT4lIAogIGJpbmRfcm93cyhtX3BhcmFtcykKCmdncGxvdChkZjEgJT4lIGZpbHRlcih0eXBlPT0ibG9zcyIpLCBhZXMoeD1hbHBoYSwgeT1iZXRhKSkgKwogIGdlb21fY29udG91cihhZXMoej1sb2cobG9zcykpLCBiaW5zID0gMTApICsKICBnZW9tX3BvaW50KGRhdGE9ZGYxICU+JSBmaWx0ZXIodHlwZT09InRyYWluaW5nIikpICsKICBnZW9tX2xpbmUoZGF0YT1kZjEgJT4lIGZpbHRlcih0eXBlPT0idHJhaW5pbmciKSkKYGBgCgpJZiB0b28gaGlnaCBhIHZhbHVlIGlzIGNob3NlbiwgdGhlIHN0ZXBzIHRha2VuIGluIGEgZ2l2ZW4gZGlyZWN0aW9uIG92ZXJzaG9vdCB0aGUgY3VydmF0dXJlLCByZXN1bHRpbmcgaW4gYW4gaW5jcmVhc2UgaW4gdGhlIGxvc3MuIElmIHRoZXNlIHN0ZXBzIGFyZSByZXBlYXRlZGx5IHRha2VuLCB0aGUgbG9zcyBkaXZlcmdlbmNlcy4KCjguIEZvciBhIHJlYXNvbmFibGUgdmFsdWUgb2YgJFxldGEkLCB2aXN1YWxpc2UgdGhlIG1lYW4gc3F1YXJlZCBsb3NzIG92ZXIgdGltZQpgYGB7cn0KbV9wYXJhbXMgPC0gZl90cmFpbl9tb2RlbCg1MDAsIDAuMDEsIHksIHgpCgpsb3NzZXMgPC0gbWFwX2RibChzZXFfYWxvbmcobV9wYXJhbXMkYWxwaGEpLCB+bXNsKHksIHgsIG1fcGFyYW1zJGFscGhhWy5dLCBtX3BhcmFtcyRiZXRhWy5dKSkKbV9wYXJhbXMgPC0gbV9wYXJhbXMgJT4lIAogIG11dGF0ZShsb3NzPWxvc3NlcykKbV9wYXJhbXMgJT4lIAogIGdncGxvdChhZXMoeD1lcG9jaCwgeT1sb3NzKSkgKwogIGdlb21fbGluZSgpCmBgYAoKOS4gT3ZlcmxheSB5b3VyIGVzdGltYXRlZCByZWdyZXNzaW9uIGxpbmUgb24gdGhlIGRhdGEKYGBge3J9Cnhfc2ltIDwtIHNlcSgtMTIsIDEyLCAwLjEpCnlfc2ltIDwtIG1fcGFyYW1zJGFscGhhW25yb3cobV9wYXJhbXMpXSArIG1fcGFyYW1zJGJldGFbbnJvdyhtX3BhcmFtcyldICogeF9zaW0KZGYgPC0gdGliYmxlKHgsIHksIHR5cGU9ImFjdHVhbCIpICU+JSAKICBiaW5kX3Jvd3ModGliYmxlKHg9eF9zaW0sIHk9eV9zaW0sIHR5cGU9InNpbXVsYXRlZCIpKQpnZ3Bsb3QoZGYgJT4lIGZpbHRlcih0eXBlPT0iYWN0dWFsIiksIGFlcyh4PXgsIHk9eSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fbGluZShkYXRhPWRmICU+JSBmaWx0ZXIodHlwZT09InNpbXVsYXRlZCIpKQoKYGBgCgoxMC4gR2VuZXJhdGUgc2ltaWxhciBkYXRhIGV4Y2VwdCB3aXRoICRcYmV0YT0wLjAyJCBhbmQgJHhcc2ltXG1hdGhjYWx7Tn0oMCwgNDAwKSQuIEhvdyBpcyB0cmFpbmluZyBpbXBhY3RlZCBmb3IgdGhpcyBtb2RlbD8gCmBgYHtyfQphbHBoYSA8LSAxCmJldGEgPC0gMC4wMgpzaWdtYSA8LSAxMApuIDwtIDEwMAp4IDwtIHJub3JtKG4sIDAsIDQwMCkKeSA8LSBhbHBoYSArIGJldGEgKiB4ICsgcm5vcm0obiwgMCwgc2lnbWEpCnBsb3QoeCwgeSkKCnBhcmFtcyA8LSBleHBhbmRfZ3JpZChhbHBoYT1zZXEoLTIsIDUsIDAuMSksIGJldGE9c2VxKC0wLjA1LCAwLjA1LCAwLjAwNSkpCnogPC0gbWFwMl9kYmwocGFyYW1zJGFscGhhLCBwYXJhbXMkYmV0YSwgfm1zbCh5LCB4LCAueCwgLnkpKQpwYXJhbXMgPC0gcGFyYW1zICU+JSAKICBtdXRhdGUobG9zcz16KQoKbV9wYXJhbXMgPC0gZl90cmFpbl9tb2RlbCgxMDAwLCAwLjAwMDAwMDEsIHksIHgpICU+JSAKICBtdXRhdGUodHlwZT0idHJhaW5pbmciKQoKZGYxIDwtIHBhcmFtcyAlPiUgCiAgbXV0YXRlKHR5cGU9Imxvc3MiKSAlPiUgCiAgYmluZF9yb3dzKG1fcGFyYW1zKQoKZ2dwbG90KGRmMSAlPiUgZmlsdGVyKHR5cGU9PSJsb3NzIiksIGFlcyh4PWFscGhhLCB5PWJldGEpKSArCiAgZ2VvbV9jb250b3VyKGFlcyh6PWxvZyhsb3NzKSksIGJpbnMgPSAxMCkgKwogIGdlb21fcG9pbnQoZGF0YT1kZjEgJT4lIGZpbHRlcih0eXBlPT0idHJhaW5pbmciKSkgKwogIGdlb21fbGluZShkYXRhPWRmMSAlPiUgZmlsdGVyKHR5cGU9PSJ0cmFpbmluZyIpKQpgYGAKClRyYWluaW5nIGlzIHZlcnkgc2xvdyBmb3IgJFxhbHBoYSQgZHVlIHRvIHRoZSBkaWZmZXJlbnQgcGFyYW1ldGVyIHNjYWxlcy4KCjExLiBTdGFuZGFyZGlzZSBib3RoIHRoZSAkeSQgYW5kICR4JCB2YXJpYWJsZXMgYW5kIHJldHVybiB0cmFpbmluZy4gSG93IGlzIHRoZSByYXRlIG9mIGNvbnZlcmdlbmNlIG5vdz8KYGBge3J9CnlfdGlsZGUgPC0gc2NhbGUoeSlbLCAxXQp4X3RpbGRlIDwtIHNjYWxlKHgpWywgMV0KCnBhcmFtcyA8LSBleHBhbmRfZ3JpZChhbHBoYT1zZXEoLTIsIDUsIDAuMSksIGJldGE9c2VxKC0xLCAyLCAwLjAxKSkKeiA8LSBtYXAyX2RibChwYXJhbXMkYWxwaGEsIHBhcmFtcyRiZXRhLCB+bXNsKHlfdGlsZGUsIHhfdGlsZGUsIC54LCAueSkpCnBhcmFtcyA8LSBwYXJhbXMgJT4lIAogIG11dGF0ZShsb3NzPXopCgptX3BhcmFtcyA8LSBmX3RyYWluX21vZGVsKDEwMDAsIDAuMDEsIHlfdGlsZGUsIHhfdGlsZGUsIGluaXRzID0gYygtMSwgLTEpKSAlPiUgCiAgbXV0YXRlKHR5cGU9InRyYWluaW5nIikKCmRmMSA8LSBwYXJhbXMgJT4lIAogIG11dGF0ZSh0eXBlPSJsb3NzIikgJT4lIAogIGJpbmRfcm93cyhtX3BhcmFtcykKCmdncGxvdChkZjEgJT4lIGZpbHRlcih0eXBlPT0ibG9zcyIpLCBhZXMoeD1hbHBoYSwgeT1iZXRhKSkgKwogIGdlb21fY29udG91cihhZXMoej1sb2cobG9zcykpLCBiaW5zID0gMTApICsKICBnZW9tX3BvaW50KGRhdGE9ZGYxICU+JSBmaWx0ZXIodHlwZT09InRyYWluaW5nIikpICsKICBnZW9tX2xpbmUoZGF0YT1kZjEgJT4lIGZpbHRlcih0eXBlPT0idHJhaW5pbmciKSkKYGBgCgpUcmFpbmluZyBub3cgbXVjaCBmYXN0ZXI6IGluIHN1bW1hcnksIGdyYWRpZW50IGRlc2NlbnQgbXVjaCBtb3JlIGVmZmljaWVudCBpZiBkYXRhIGFyZSBvbiBzaW1pbGFyIHNjYWxlcy4gU28sIGZlYXR1cmUgc2NhbGluZyByZWFsbHkgaW1wb3J0YW50LgoKMTIuIEdlbmVyYXRlIHJlZ3Jlc3Npb24gZGF0YSB2aWE6CgpcYmVnaW57ZXF1YXRpb259CnlfaSBcc3RhY2tyZWx7aS5pLmQufXtcc2ltfSBcbWF0aGNhbHtOfShcYWxwaGEgKyBcYmV0YSB4X2ksIFxzaWdtYSkKXGVuZHtlcXVhdGlvbn0KCndoZXJlICRcYWxwaGE9MSQsICRcYmV0YT0xJCwgJFxzaWdtYT01JCBhbmQgJHhfaSBcc3RhY2tyZWx7aS5pLmQufXtcc2ltfVxtYXRoY2Fse059KDAsIDQpJCBmb3IgJGkgPTEsLi4uLDMwJC4KClVzZSBncmFkaWVudCBkZXNjZW50IHRvIG9idGFpbiBwYXJhbWV0ZXIgZXN0aW1hdGVzIGZvciBhIHF1YWRyYXRpYyByZWdyZXNzaW9uIGFuZCBncmFwaCB0aGUgZml0dGVkIGxpbmUgdmVyc3VzIHRoZSBkYXRhLgoKYGBge3J9CmFscGhhIDwtIDEKYmV0YSA8LSAxCnNpZ21hIDwtIDUKbiA8LSAzMAp4IDwtIHJub3JtKG4sIDAsIDQpCnkgPC0gYWxwaGEgKyBiZXRhICogeCArIHJub3JtKG4sIDAsIHNpZ21hKQoKdXBkYXRlX3BhcmFtZXRlcnNfcXVhZCA8LSBmdW5jdGlvbihhbHBoYSwgYmV0YSwgZ2FtbWEsIHksIHgsIGV0YSkgewogIEsgPC0gbGVuZ3RoKHkpCiAgZGxfZGEgPC0gLSgyIC8gSykgKiBzdW0oeSAtIChhbHBoYSArIGJldGEgKiB4ICsgZ2FtbWEgKiB4XjIpKQogIGRsX2RiIDwtIC0oMiAvIEspICogc3VtKHggKiAoeSAtIChhbHBoYSArIGJldGEgKiB4ICsgZ2FtbWEgKiB4XjIpKSkKICBkbF9kZyA8LSAtKDIgLyBLKSAqIHN1bSh4XjIgKiAoeSAtIChhbHBoYSArIGJldGEgKiB4ICsgZ2FtbWEgKiB4XjIpKSkKICAKICBhbHBoYSA8LSBhbHBoYSAtIGV0YSAqIGRsX2RhCiAgYmV0YSA8LSBiZXRhIC0gZXRhICogZGxfZGIKICBnYW1tYSA8LSBnYW1tYSAtIGV0YSAqIGRsX2RnCiAgcmV0dXJuKGMoYWxwaGEsIGJldGEsIGdhbW1hKSkKfQoKbl9lcG9jaHMgPC0gNTAwMAptX3BhcmFtcyA8LSBtYXRyaXgobnJvdyA9IG5fZXBvY2hzLCBuY29sID0gMykKbV9wYXJhbXNbMSwgXSA8LSBjKDAsIDAsIDApCmZvcihpIGluIDI6bl9lcG9jaHMpIHsKICBtX3BhcmFtc1tpLCBdIDwtIHVwZGF0ZV9wYXJhbWV0ZXJzX3F1YWQobV9wYXJhbXNbKGkgLSAxKSwgMV0sIG1fcGFyYW1zWyhpIC0gMSksIDJdLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtX3BhcmFtc1soaSAtIDEpLCAzXSwgeSwgeCwgMC4wMDEpCn0KbV9wYXJhbXMgPC0gYXMuZGF0YS5mcmFtZShtX3BhcmFtcykKY29sbmFtZXMobV9wYXJhbXMpIDwtIGMoImFscGhhIiwgImJldGEiLCAiZ2FtbWEiKQptX3BhcmFtcyA8LSBtX3BhcmFtcyAlPiUgCiAgbXV0YXRlKGVwb2NoPXNlcV9hbG9uZyhhbHBoYSkpCm1fcGFyYW1zICU+JSAKICBtZWx0KGlkLnZhcnM9ImVwb2NoIikgJT4lIAogIGdncGxvdChhZXMoeD1lcG9jaCwgeT12YWx1ZSwgY29sb3VyPWFzLmZhY3Rvcih2YXJpYWJsZSkpKSArCiAgZ2VvbV9saW5lKCkKCnhfc2ltIDwtIHNlcSgtMTIsIDEyLCAwLjEpCnlfc2ltIDwtIG1fcGFyYW1zJGFscGhhW25yb3cobV9wYXJhbXMpXSArIG1fcGFyYW1zJGJldGFbbnJvdyhtX3BhcmFtcyldICogeF9zaW0gKyBtX3BhcmFtcyRnYW1tYVtucm93KG1fcGFyYW1zKV0gKiB4X3NpbV4yCgpmX3F1YWRyYXRpY19yZWdyZXNzaW9uIDwtIGZ1bmN0aW9uKHgsIGFscGhhLCBiZXRhLCBnYW1tYSkgewogIHJldHVybihhbHBoYSArIGJldGEgKiB4ICsgZ2FtbWEgKiB4XjIpCn0KYSA8LSBtX3BhcmFtcyRhbHBoYVtucm93KG1fcGFyYW1zKV0KYiA8LSBtX3BhcmFtcyRiZXRhW25yb3cobV9wYXJhbXMpXQpjIDwtIG1fcGFyYW1zJGdhbW1hW25yb3cobV9wYXJhbXMpXQpmX3F1YWQgPC0gZnVuY3Rpb24oeCkgZl9xdWFkcmF0aWNfcmVncmVzc2lvbih4LCBhLCBiLCBjKQoKZGYgPC0gdGliYmxlKHgsIHksIHR5cGU9ImFjdHVhbCIpICU+JSAKICBiaW5kX3Jvd3ModGliYmxlKHg9eF9zaW0sIHk9eV9zaW0sIHR5cGU9InNpbXVsYXRlZCIpKQpnZ3Bsb3QoZGYgJT4lIGZpbHRlcih0eXBlPT0iYWN0dWFsIiksIGFlcyh4PXgsIHk9eSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fbGluZShkYXRhPWRmICU+JSBmaWx0ZXIodHlwZT09InNpbXVsYXRlZCIpKQpgYGAKCjEzLiBGaXQgdGhlIHNhbWUgZGF0YSBhcyBnZW5lcmF0ZWQgaW4gdGhlIHByZXZpb3VzIHF1ZXN0aW9uIHVzaW5nIGxpbmVhciByZWdyZXNzaW9uLiBDb21wYXJlIHRoZSBmaXQgb2YgdGhlIHF1YWRyYXRpYyBhbmQgbGluZWFyIG1vZGVscyBvbiBuZXcgZGF0YSBnZW5lcmF0ZWQgZnJvbSB0aGUgc2FtZSBkYXRhIGdlbmVyYXRpbmcgcHJvY2VzcyBleGNlcHQ6CgokeF9pIFxzdGFja3JlbHtpLmkuZC59e1xzaW19XG1hdGhjYWx7Tn0oNSwgNCkkIGZvciAkaSA9MSwuLi4sMzAkCgpXaGljaCBoYXMgYSBiZXR0ZXIgcm9vdCBtZWFuIHNxdWFyZWQgYWNjdXJhY3kgb24gdGhlIG5ldyBkYXRhPyBXaHk/CgpgYGB7cn0KbV9wYXJhbXMgPC0gZl90cmFpbl9tb2RlbCgyMDAwLCAwLjAxLCB5LCB4KQoKZl9saW5lYXJfcmVncmVzc2lvbiA8LSBmdW5jdGlvbih4LCBhbHBoYSwgYmV0YSkgewogIHJldHVybihhbHBoYSArIGJldGEgKiB4KQp9CgpmX2xpbmVhciA8LSBmdW5jdGlvbih4KSBmX2xpbmVhcl9yZWdyZXNzaW9uKHgsbV9wYXJhbXMkYWxwaGFbbnJvdyhtX3BhcmFtcyldLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtX3BhcmFtcyRiZXRhW25yb3cobV9wYXJhbXMpXSkKCngxIDwtIHJub3JtKG4sIDUsIDQpCnkxIDwtIGFscGhhICsgYmV0YSAqIHgxICsgcm5vcm0obiwgMCwgc2lnbWEpCgojIHByZWRpY3QgeSBhbmQgY2FsY3VsYXRlIFJNU0VzCnkxX2xpbmVhciA8LSBtYXBfZGJsKHgxLCBmX2xpbmVhcikKeTFfcXVhZCA8LSBtYXBfZGJsKHgxLCBmX3F1YWQpCgpzcXJ0KG1lYW4oeTFfbGluZWFyIC0geTEpXjIpCnNxcnQobWVhbih5MV9xdWFkIC0geTEpXjIpCmBgYAoKVGhlIGxpbmVhciBtb2RlbCBnZW5lcmFsaXNlcyBiZXR0ZXIuIFRoZSBxdWFkcmF0aWMgbW9kZWwgaGFzIG92ZXJmaXR0ZWQgdGhlIHRyYWluaW5nIGRhdGEsIHNvIGRvZXNuJ3QgZXh0cmFwb2xhdGUgd2VsbCB0byBuZXcgZGF0YSByZWdpbWVzLgoKIyBLTk4gcmVncmVzc2lvbgoKMS4gR2VuZXJhdGUgJG49MjAwJCBwYWlyZWQgJCh4X2ksIHlfaSkkIGRhdGEgcG9pbnRzIGJ5IHNhbXBsaW5nOgoKXGJlZ2lue2VxdWF0aW9ufQp4X2kgXHNpbSBcbWF0aGNhbHtOfSgwLCA0KQpcZW5ke2VxdWF0aW9ufQoKYW5kIHRoZW46CgpcYmVnaW57ZXF1YXRpb259CnlfaSBcc2ltIFxtYXRoY2Fse059KHNpbih4X2kpLCAwLjIpClxlbmR7ZXF1YXRpb259CgpgYGB7cn0KbGlicmFyeShSQU5OKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbiA8LSAyMDAKeCA8LSBybm9ybShuLCAwLCA0KQp5IDwtIHNpbih4KSArIHJub3JtKG4sIDAsIDAuMikKCmRmIDwtIHRpYmJsZSh4LCB5KSAlPiUgCiAgbXV0YXRlKHR5cGU9ImFjdHVhbCIpCgpnZ3Bsb3QoZGYgJT4lIGZpbHRlcih0eXBlPT0iYWN0dWFsIiksIGFlcyh4PXgsIHk9eSkpICsKICBnZW9tX3BvaW50KCkKYGBgCgoyLiBDcmVhdGUgYSBmdW5jdGlvbiB3aGljaCByZXR1cm5zIHRoZSBLTk5zIGZvciBhIGdpdmVuICR4JCB2YWx1ZSAoYW5kICRrJCkKYGBge3J9CiMgdXNlIFJBTk4gbGlicmFyeSAod2hpY2ggaW1wbGVtZW50cyBLRCB0cmVlcykKeF90aWxkZSA8LSAtMi41CmZpdCA8LSBubjIoeCwgYyh4X3RpbGRlKSwgayA9IDEwKQppZHhzIDwtIGZpdCRubi5pZHhbMSwgXQoKIyB2aXN1YWxpc2UgbmVhcmVzdCBuZWlnaGJvdXJzIGFyb3VuZCBhIHBvaW50CnRpYmJsZSh4LCB5KSAlPiUgCiAgbXV0YXRlKGlkeD1zZXFfYWxvbmcoeCkpICU+JSAKICBtdXRhdGUobm49aWR4JWluJWlkeHMpICU+JSAKICBnZ3Bsb3QoYWVzKHg9eCwgeT15LCBjb2xvdXI9YXMuZmFjdG9yKG5uKSkpICsKICBnZW9tX3BvaW50KCkKYGBgCgozLiBVc2luZyB5b3VyIGFuc3dlciB0byB0aGUgcHJldmlvdXMgY29kZSwgY3JlYXRlIGEgS05OIHJlZ3Jlc3Npb24gZXN0aW1hdG9yIHVzaW5nIEV1Y2xpZGVhbiBkaXN0YW5jZSB3aXRoIHZhcmlvdXMgdmFsdWVzIG9mICRrJCB0byBwcmVkaWN0ICR5JC4gSG93IGRvZXMgJGskIGluZmx1ZW5jZSB0aGUgcmVzdWx0cz8KYGBge3J9CmZfa2tfcmVncmVzc2lvbiA8LSBmdW5jdGlvbih4X3RpbGRlLCB4LCB5LCBrPTEwKSB7CiAgZml0IDwtIG5uMih4LCBjKHhfdGlsZGUpLCBrID0gaykKICBpZHhzIDwtIGZpdCRubi5pZHhbMSwgXQogIHJldHVybihtZWFuKHlbaWR4c10pKQp9Cgp4X3NpbSA8LSBzZXEoLTEwLCAxMCwgMC4xKQoKIyBrID0gMQp5X3NpbSA8LSBtYXBfZGJsKHhfc2ltLCB+Zl9ra19yZWdyZXNzaW9uKC4sIHgsIHksIGs9MTApKQoKZGYgPC0gdGliYmxlKHgsIHkpICU+JSAKICBtdXRhdGUodHlwZT0iYWN0dWFsIikgJT4lIAogIGJpbmRfcm93cyh0aWJibGUoeD14X3NpbSwgeT15X3NpbSkgJT4lIG11dGF0ZSh0eXBlPSJyZWdyZXNzaW9uIikpCgpnZ3Bsb3QoZGYgJT4lIGZpbHRlcih0eXBlPT0iYWN0dWFsIiksIGFlcyh4PXgsIHk9eSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fbGluZShkYXRhPWRmICU+JSBmaWx0ZXIodHlwZT09InJlZ3Jlc3Npb24iKSwgY29sb3VyPSJibHVlIikKCiMgayA9IDEKeV9zaW0gPC0gbWFwX2RibCh4X3NpbSwgfmZfa2tfcmVncmVzc2lvbiguLCB4LCB5LCBrPTEpKQoKZGYgPC0gdGliYmxlKHgsIHkpICU+JSAKICBtdXRhdGUodHlwZT0iYWN0dWFsIikgJT4lIAogIGJpbmRfcm93cyh0aWJibGUoeD14X3NpbSwgeT15X3NpbSkgJT4lIG11dGF0ZSh0eXBlPSJyZWdyZXNzaW9uIikpCgpnZ3Bsb3QoZGYgJT4lIGZpbHRlcih0eXBlPT0iYWN0dWFsIiksIGFlcyh4PXgsIHk9eSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fbGluZShkYXRhPWRmICU+JSBmaWx0ZXIodHlwZT09InJlZ3Jlc3Npb24iKSwgY29sb3VyPSJibHVlIikKCiMgayA9IDEwMAp5X3NpbSA8LSBtYXBfZGJsKHhfc2ltLCB+Zl9ra19yZWdyZXNzaW9uKC4sIHgsIHksIGs9MTAwKSkKCmRmIDwtIHRpYmJsZSh4LCB5KSAlPiUgCiAgbXV0YXRlKHR5cGU9ImFjdHVhbCIpICU+JSAKICBiaW5kX3Jvd3ModGliYmxlKHg9eF9zaW0sIHk9eV9zaW0pICU+JSBtdXRhdGUodHlwZT0icmVncmVzc2lvbiIpKQoKZ2dwbG90KGRmICU+JSBmaWx0ZXIodHlwZT09ImFjdHVhbCIpLCBhZXMoeD14LCB5PXkpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2xpbmUoZGF0YT1kZiAlPiUgZmlsdGVyKHR5cGU9PSJyZWdyZXNzaW9uIiksIGNvbG91cj0iYmx1ZSIpCmBgYAoK